T0 Summer Stress ANOVAs

Output contents

• Final figures & summary stats for Part I pub, ANOVA models:
• ANOVA boxplots/means (Pooled Site, Topography, Microhabitat Type) - pooled across all levels.
  
• One-way ANOVA families (Stress & Performance ~ Site, ~ Topo, ~ Microhabitat):
  • Histograms - normality assessment raw data
    
  • Homogeneity Tests raw data
  • Anovas & residual plots
  • Nonparametric (Kruskal-Wallis) if needed.
  • Family p-value adjustments (or not? I am against these in principle, but especially for exploratory analysis like this ~ no controlled experiments here, so messy data is the norm!)

Notes about this data

• Data was collected by T. Clark & A. Russel in the Sheep Range (DNWR), Nevada.
  
• Experiment = Summer stress Assay on n=96 microsites (minus 2 wrong species) = n=94 sample units.
• Sample units = high-density Syntrichia caninervis biocrust microsites (20 x 20 cm).
• Time = Time zero (T0) at which point physiological measurements were taken (30min post-rehydration, before significant desiccation recovery had occurred).
• Script error correction for Phi-PSII (PhiPSII): generic correction applied to Assay 1 samples to fix light-adapted photosynthesis measurements accidentily made at lower light (PAR = 10) instead of PAR = 150. See below for details. Corrected variable is called PhiPSII.
• Sunseeker shade index is included in this dataset but will be removed before upload to GitHub. Solarveg is used instead (on continuous models in another file).

Set Gobal Chunk Options & Load libraries

• WARNING: Setting cache = TRUE can cause random errors in the code (like variables that cannot be found suddenly). If error arises, set cache=FALSE and rerun.

knitr::opts_chunk$set(
    echo = T,
    warning = FALSE,
    cache = FALSE,
    collapse = T
)
getwd()
library(tidyverse) #WHY WAS THIS NOT LOADED PREVIOUSLY??!!!
#library(nlme) #mixed modeling - delete if not used
library(plyr)#revalue() factors and reorder them! dplyr::recode not working
library(dplyr)  #summarise_at()
library(ggplot2) #graphing
#library(lattice) #delete if not used
#library(gtable)
library(grid) #tables?
library(gridExtra)#masks combine() from dplyr
library(pander) #table output
library(Publish) #delete if not used
library(ggpubr) #combine plots
library(kableExtra)
library(devtools) #
library(broom) #
library(purrr)#
#Modeling libraries
library(car)#sophisticated anova() modeling
#library(rstatix)#anova_test() - masks mutate() and desc() from plyr package!
library(sjstats)#anova_stats()
library(pwr)#unknown function but the sjstats package must use it to calculate Cohen's D!
#library(MDmisc)#download from GitHub - 

Import Data & Set Factor Levels

• Nooooooooo the era of automatic factors is over.
• Question: How can you convert a bunch of columns to factors simultaneously?

#rm(Recovery)
D<-read.csv("T0_Stress_R.csv")

# Set factors (no longer automatic)
D$Moss<-as.factor(D$Moss)
D$Microsite<-as.factor(D$Microsite)
D$Habitat<-as.factor(D$Habitat)

#Remove other species & drop PAR rejected for tissue sampling
D<-droplevels(D[!c(D$Microsite=="PJ-N-91" | D$Microsite =="PJ-S-67a" | 
                    D$Microsite=="BB-N-35" | D$Microsite=="BB-S-53"),])
length(levels(D$Microsite))#Down to 94 units

# Order factor levels
D$Topography<-factor(D$Topography, levels=c("S", "F", "N"))
D$Veg.Type<-factor(D$Veg.Type, levels=c("CB", "BB", "PJ"))
D$Veg.Type<-plyr::revalue(D$Veg.Type, c("CB" = "Low-scrubland", "BB" = "Mid-shrubland", "PJ" ="High-woodland"))

D$Site<- factor(D$Site, levels=c("1","2","3"))
#levels(D$Site)<- list("1" = "Low-scrubland", "2"= "Mid-shrubland", "3" = "High-woodland")
## [1] 94

Script error correction (PhiPSII, qP)

• Opted to not do a correction, but just deleted Assay 1 for PhiPSII analysis.
• Tested the correlation between Fv/Fm and Phi-PSII for Assays 2, 3, & 4. Select a correction coefficient for Assay 1 that creates this level of correlation for Assay 1.
• Recommended correction of -0.22 per light curves created several (five?) negative values, so I choose the closest correction coefficient that kept values positive = -0.17.
• Question: why is Assay 1 variance SO much higher than the other assays?? Nick project…
• Question: what is expected correlation between Fv/Fm and Phi-PSII? Assay 2 - 4 were about 0.91, while assay 1 cor = 0.96, so too high? (Correlation doesn’t change after the linear correction.)

#Mean PhiPSII (Phi-PSII) of the Assays (factor levels)
D %>% group_by(., Assay) %>% summarise_at(vars(PhiPSII), mean)

#Mean difference (Assay 1 - the others)
with(D, mean(PhiPSII[Assay!=1]))#0.2464286
## [1] 0.2464286
with(D, mean(PhiPSII[Assay!=1])) - with(D, mean(PhiPSII[Assay==1]))#-0.1834881 
## [1] -0.1834881
#D<-D %>% mutate(PhiPSII = ifelse(Assay ==1, PhiPSII-0.17, PhiPSII))# Correction = -0.183
#D$Moss[D$PhiPSII< 0]#Chosen correction where no values become negative in Assay 1.
#D<-D %>% mutate(qPhoto = ifelse(Assay ==1, qP-0.22, qP))
#D$Moss[D$qPhoto<0]#Worked all are positive

#Correlation of Fv/Fm for Assay 2 - 4
#Assay 2 = 0.9056776, Assay 3 = 0.9093973, Assay 4 = 0.9148743, Assay 1 = 0.9627869
D %>% group_by(., Assay) %>% do(as.data.frame(cor(.$Fv.Fm, .$PhiPSII), data=.))

#Correlation of the pooled dataset, thus, improves after correction, but too much?
with(D, cor(PhiPSII, Fv.Fm))#0.74
## [1] 0.7444226
#with(D, cor(PhiPSII, Fv.Fm))#0.909 

Delete Assay 1

• How to remove certain rows for subsequent analysis?
• Maybe just opt to remove them during each analysis or plot using arguments to drop certain rows?

#See change after correction medians align!
par(mar = c(4, 4, .1, .1))
with(D, boxplot(PhiPSII~Assay))
with(D[D$Assay != 1,], boxplot(PhiPSII~Assay))#Dropped Assay 1
# Density plots;transparent
ggplot(D[D$Assay != 1,], aes(x=PhiPSII, fill=Veg.Type)) + geom_density(alpha=.6)#

#————————- # SUMMARY STATS ## Function * New way to enter functions after 2019 packages update:

summary(D)# means, min, max, etc for all variables
##       Time        Moss     Moss.number      Moss.sub            Microsite 
##  Min.   :0   10     : 1   Min.   : 1.00   Length:94          BB-F-37 : 1  
##  1st Qu.:0   11     : 1   1st Qu.:22.25   Class :character   BB-F-38a: 1  
##  Median :0   12     : 1   Median :47.50   Mode  :character   BB-F-39a: 1  
##  Mean   :0   13a    : 1   Mean   :47.12                      BB-F-40a: 1  
##  3rd Qu.:0   14a    : 1   3rd Qu.:71.75                      BB-F-41a: 1  
##  Max.   :0   15a    : 1   Max.   :96.00                      BB-F-42a: 1  
##              (Other):88                                      (Other) :88  
##  Old.Habitat              Habitat       Index       Site            Veg.Type 
##  Length:94          Canopy    :48   Min.   : 3.00   1:26   Low-scrubland:26  
##  Class :character   Dripline  :28   1st Qu.: 8.00   2:34   Mid-shrubland:34  
##  Mode  :character   Interspace:18   Median :14.00   3:34   High-woodland:34  
##                                     Mean   :14.44                            
##                                     3rd Qu.:20.00                            
##                                     Max.   :26.00                            
##                                                                              
##  Topography  Shade.Index         PDIR            Slope            Aspect      
##  S:23       Min.   : 6.00   Min.   :0.7717   Min.   : 0.000   Min.   :  0.00  
##  F:37       1st Qu.:12.00   1st Qu.:0.9142   1st Qu.: 2.000   1st Qu.: 90.75  
##  N:34       Median :20.00   Median :0.9846   Median : 7.000   Median :200.00  
##             Mean   :17.87   Mean   :0.9486   Mean   : 7.468   Mean   :194.88  
##             3rd Qu.:23.00   3rd Qu.:0.9997   3rd Qu.:12.000   3rd Qu.:270.00  
##             Max.   :27.00   Max.   :1.0496   Max.   :18.000   Max.   :357.00  
##                                                               NA's   :4       
##    VegPercent        SolarVeg          Camera           Pair     
##  Min.   :  0.00   Min.   :  0.00   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.: 22.22   1st Qu.: 29.17   1st Qu.: 6.25   1st Qu.: 9.5  
##  Median : 38.89   Median : 43.75   Median :11.50   Median :19.0  
##  Mean   : 40.66   Mean   : 45.26   Mean   :11.73   Mean   :18.9  
##  3rd Qu.: 60.41   3rd Qu.: 62.50   3rd Qu.:16.75   3rd Qu.:28.5  
##  Max.   :100.00   Max.   :100.00   Max.   :23.00   Max.   :37.0  
##                                    NA's   :72      NA's   :55    
##  Notes.Urgent        Assay.Set             Assay            Clip      
##  Length:94          Length:94          Min.   :1.000   Min.   : 1.00  
##  Class :character   Class :character   1st Qu.:1.250   1st Qu.: 6.25  
##  Mode  :character   Mode  :character   Median :2.000   Median :12.00  
##                                        Mean   :2.468   Mean   :12.32  
##                                        3rd Qu.:3.000   3rd Qu.:18.00  
##                                        Max.   :4.000   Max.   :24.00  
##                                                                       
##        Fo               Fm              Fv             Fv.Fm       
##  Min.   : 30.00   Min.   : 52.0   Min.   : 12.00   Min.   :0.1300  
##  1st Qu.: 79.25   1st Qu.:129.5   1st Qu.: 45.25   1st Qu.:0.3515  
##  Median : 96.00   Median :183.5   Median : 93.50   Median :0.4770  
##  Mean   : 99.02   Mean   :195.0   Mean   : 95.97   Mean   :0.4546  
##  3rd Qu.:118.50   3rd Qu.:248.8   3rd Qu.:132.00   3rd Qu.:0.5697  
##  Max.   :174.00   Max.   :419.0   Max.   :265.00   Max.   :0.7370  
##                                                                    
##        Fs              Fm.L            Fo.L             Fv.L       
##  Min.   : 29.00   Min.   : 46.0   Min.   : 17.00   Min.   : 15.00  
##  1st Qu.: 76.25   1st Qu.:101.0   1st Qu.: 52.50   1st Qu.: 42.25  
##  Median : 92.50   Median :137.5   Median : 71.50   Median : 61.00  
##  Mean   : 96.63   Mean   :141.8   Mean   : 70.89   Mean   : 70.86  
##  3rd Qu.:118.00   3rd Qu.:175.5   3rd Qu.: 89.25   3rd Qu.: 95.00  
##  Max.   :190.00   Max.   :319.0   Max.   :149.00   Max.   :211.00  
##                                                                    
##     Fv.Fm.L          PhiPSII             qP              qNP         
##  Min.   :0.2230   Min.   :0.0430   Min.   :0.1330   Min.   :-0.8570  
##  1st Qu.:0.4050   1st Qu.:0.1900   1st Qu.:0.4903   1st Qu.: 0.1765  
##  Median :0.4905   Median :0.2685   Median :0.5820   Median : 0.4100  
##  Mean   :0.4881   Mean   :0.2933   Mean   :0.5826   Mean   : 0.3651  
##  3rd Qu.:0.5780   3rd Qu.:0.3635   3rd Qu.:0.6707   3rd Qu.: 0.5870  
##  Max.   :0.7360   Max.   :0.6930   Max.   :0.9460   Max.   : 0.7970  
##                                                                      
##       NPQ             Shrubindex        Topoindex        Siteindex       
##  Min.   :-0.19100   Min.   :-1.0000   Min.   :-1.000   Min.   :-1.00000  
##  1st Qu.: 0.08575   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:-1.00000  
##  Median : 0.32000   Median : 1.0000   Median : 0.000   Median : 0.00000  
##  Mean   : 0.43293   Mean   : 0.3191   Mean   : 0.117   Mean   : 0.08511  
##  3rd Qu.: 0.68800   3rd Qu.: 1.0000   3rd Qu.: 1.000   3rd Qu.: 1.00000  
##  Max.   : 2.54200   Max.   : 1.0000   Max.   : 1.000   Max.   : 1.00000  
##                                                                          
##   Multibuffer     
##  Min.   :-2.0000  
##  1st Qu.: 0.0000  
##  Median : 1.0000  
##  Mean   : 0.5213  
##  3rd Qu.: 1.0000  
##  Max.   : 3.0000  
## 

#Stats needed for plotting
SE<- function(x) round((sqrt(var(x)/length(x))),3)
Myfunctions = list (
  Mean = function(x) round(mean(x),3),
  Max= function(x)round(max(x),3),
  Min= function(x)round(min(x),3),
  Range= function(x) round((max(x)-min(x)),3),
  SE= function(x) round((sqrt(var(x)/length(x))),3),
  Upper= function(x) round(mean(x)+SE(x),3),
  Lower= function(x) round(mean(x)-SE(x),3),
  Var = function(x) round(var(x),3),
  SD = function(x) round(sd(x),3)
)

Tables

###Stress (Fv/Fm) & Performance (PhiPSII)

#Fv/Fm Tables
#Site table
T.Site.Fvfm<-D %>% 
  group_by(Site) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Site.Fvfm, caption = "Fv/Fm Stats by Site", format="pandoc", padding=2)

Fv/Fm Stats by Site

Site	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
1	0.440	0.602	0.183	0.419	0.021	0.461	0.419	0.011	0.105
2	0.429	0.737	0.130	0.607	0.027	0.456	0.402	0.026	0.160
3	0.491	0.686	0.144	0.542	0.028	0.519	0.463	0.026	0.163

#Topography tables
T.Topo.Fvfm<-D %>% 
  group_by(Topography) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Topo.Fvfm, caption = "Fv/Fm Stats by Topography", format="pandoc", padding=2)

Fv/Fm Stats by Topography

Topography	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
S	0.432	0.656	0.20	0.456	0.026	0.458	0.406	0.016	0.125
F	0.391	0.647	0.13	0.517	0.027	0.418	0.364	0.026	0.162
N	0.539	0.737	0.25	0.487	0.018	0.557	0.521	0.012	0.107

#Microhabitat tables
T.Hab.Fvfm<-D %>% 
  group_by(Habitat) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Fv/Fm Stats by Habitat", format="pandoc", padding=2)

Fv/Fm Stats by Habitat

Habitat	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
Canopy	0.516	0.737	0.174	0.563	0.020	0.536	0.496	0.018	0.136
Dripline	0.421	0.647	0.130	0.517	0.028	0.449	0.393	0.021	0.146
Interspace	0.342	0.508	0.160	0.348	0.025	0.367	0.317	0.011	0.105

#Phi-PSII Tables
##Site Table
T.Site.Phi<-D %>% 
  group_by(Site) %>% 
  summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Site.Phi, caption = "PHi-PSII Stats by Site", format="pandoc", padding=2)

PHi-PSII Stats by Site

Site	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
1	0.256	0.520	0.150	0.37	0.019	0.275	0.237	0.010	0.098
2	0.281	0.693	0.043	0.65	0.024	0.305	0.257	0.020	0.141
3	0.334	0.622	0.072	0.55	0.026	0.360	0.308	0.023	0.151

#Topo table
T.Topo.Phi<-D %>% 
  group_by(Topography) %>% 
  summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Topo.Phi, caption = "PHi-PSII Stats by Topography", format="pandoc", padding=2)

PHi-PSII Stats by Topography

Topography	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
S	0.301	0.585	0.098	0.487	0.026	0.327	0.275	0.016	0.127
F	0.238	0.621	0.043	0.578	0.020	0.258	0.218	0.015	0.122
N	0.348	0.693	0.072	0.621	0.024	0.372	0.324	0.019	0.139

#Microhabitat table
T.Hab.Phi<-D %>% 
  group_by(Habitat) %>% 
  summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Phi-PSII Stats by Habitat", format="pandoc", padding=2)

Phi-PSII Stats by Habitat

Habitat	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
Canopy	0.516	0.737	0.174	0.563	0.020	0.536	0.496	0.018	0.136
Dripline	0.421	0.647	0.130	0.517	0.028	0.449	0.393	0.021	0.146
Interspace	0.342	0.508	0.160	0.348	0.025	0.367	0.317	0.011	0.105
#————-	———	–
# BOXPLOTS & M	EANS

Life zone site boxplots

• Outlier (Phi-PSII) is Moss 34 at Site 2. How could this reading be SO high? Was Fv/Fm high also? Yes: 0.737.
• Question: what dimensions to save for publication? (fig.width=4, fig.height = 1.7)

#Fv/Fm life zone boxplots
PlotA<-ggplot(D, aes(x=(Site), y=Fv.Fm))+
  geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
  geom_point(aes(fill=(Site)), position = position_dodge2(width=0.3), size=1, 
             color="black")+
  #stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
  geom_point(data = T.Site.Fvfm, aes(x=Site, y=Mean), color="white")+
  geom_errorbar(data=T.Site.Fvfm, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
  theme_linedraw()+
  labs(y='\u2190 Stress (Fv/Fm)', x="Elevation life-zone (site)")+
  #scale_x_discrete(labels=c("1" = "Low-scrubland", "2" = "Mid-shrubland", "3"= "High-woodland"))+
  scale_x_discrete(labels=c("1" = "Low-site", "2" = "Mid-site", "3"= "High-site"))+
  scale_y_continuous(breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75))+
  theme(legend.position = "none", text=element_text(size=15), axis.title.x = 
          element_blank(), panel.grid.minor = element_blank())+
  scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))# panel.grid.major = element_blank(),
#PlotA

#Phi-PSII life zone boxplots
PlotB<-ggplot(D, aes(x=(Site), y=PhiPSII))+
  geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
  geom_point(aes(fill=(Site)), position = position_dodge2( width=0.3), size=1, 
             color="black")+
  #stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
  geom_point(data = T.Site.Phi, aes(x=Site, y=Mean), color="white")+
  geom_errorbar(data=T.Site.Phi, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
  theme_linedraw()+
  labs(y='Performance (\u03A6PSII)', x="Elevation life-zone (site)")+
  scale_x_discrete(labels=c("1" = "Low-site", "2" = "Mid-site", "3"= "High-site"))+
  scale_y_continuous(breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75))+
  theme(legend.position = "none", text=element_text(size=15), axis.title.x = 
          element_blank(), panel.grid.minor = element_blank())+
  scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))# panel.grid.major = element_blank(),
#PlotB

#Plot together
PlotsAB <-ggarrange(PlotA, PlotB, labels=c("a", "b"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsAB<- annotate_figure(PlotsAB, top = text_grob("Summer stress assay vs. elevation-life zone sites", size = 15, rot=0, just= "center")) #Add titles
PlotsAB

Topography Boxplots

#Fv/Fm topo boxplots
PlotC<-ggplot(D)+
    aes(x=Topography, y = Fv.Fm, fill = Topography) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Topo.Fvfm, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper), 
                width=.2, size=0.7, color="white") +
    scale_y_continuous(name = "\u2190 Stress signal (Fv/Fm)",breaks = seq(0, 0.7, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Topography (pooled life zones)", 
                     labels=c("S-facing", "Flat", "N-facing")) +
   #ggtitle("Summer stress assay (Series 1.1)") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
       theme(panel.background = element_rect(fill = "white"), axis.title.x = element_blank(), panel.grid.minor = element_blank())
#PlotC

#Phi topo boxplots
PlotD<-ggplot(D)+
    aes(x=Topography, y = PhiPSII, fill = Topography) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Topo.Phi, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper), 
                width=.2, size=0.7, color="white") +
    scale_y_continuous(name = "Performance (\u03A6PSII)",breaks = seq(0, 0.7, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Topography (pooled life zones)", 
                     labels=c("S-facing", "Flat", "N-facing")) +
   #ggtitle("Summer stress assay (Series 1.1)") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
       theme(panel.background = element_rect(fill = "white"), axis.title.x = element_blank(), panel.grid.minor = element_blank())
#PlotD

#Plot together
PlotsCD <-ggarrange(PlotC, PlotD, labels=c("c", "d"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsCD<- annotate_figure(PlotsCD, top = text_grob("Summer stress assay vs. Topography zones", size = 15, rot=0, just= "center")) #Add titles
PlotsCD

Microhabitat boxplots

#Fv/fm
PlotE<-ggplot(D)+
    aes(x=Habitat, y =Fv.Fm, fill = Habitat) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Hab.Fvfm, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper), 
                width=.15, size=0.7, color= "white") +
    scale_y_continuous(name = "\u2190 Stress signal (Fv/Fm)",breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Microhabitat type") +
  # ggtitle("Micro-habitat buffering") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("darkgreen", "yellowgreen", "tan"))+
       theme(panel.background = element_rect(fill = "white"), panel.grid.minor = element_blank(),axis.title.x = element_blank())
#PlotE 

#Phi
PlotF<-ggplot(D)+
    aes(x=Habitat, y =PhiPSII, fill = Habitat) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Hab.Phi, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper), 
                width=.15, size=0.7, color= "white") +
    scale_y_continuous(name = "Performance (\u03A6PSII)", breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Microhabitat type") +
  # ggtitle("Micro-habitat buffering") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("darkgreen", "yellowgreen", "tan"))+
       theme(panel.background = element_rect(fill = "white"), panel.grid.minor = element_blank(),axis.title.x = element_blank())
#PlotF 

#Plot together
PlotsEF <-ggarrange(PlotE, PlotF, labels=c("e", "f"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsEF<- annotate_figure(PlotsEF, top = text_grob("Summer stress assay vs. Habitat types", size = 15, rot=0, just= "center")) #Add titles
PlotsEF

Final boxplots: Stress & performance vs 3 scales of buffering

#Plot A-F together
PlotsA_F <-ggarrange(PlotA, PlotB, PlotC, PlotD, PlotE, PlotF, labels=c("a","b","c","d","e", "f"), widths=1:1, ncol=2, nrow=3, font.label= list(size=20)) #Add letters & up the size
PlotsA_F<- annotate_figure(PlotsA_F, top = text_grob("Summer stress & performance vs. 3 scales of habitat buffering", size = 16, rot=0, just= "center")) #Add bottom/top titles
PlotsA_F

#——————— # VARIANCE TESTS

• Goal: Test null of equal variance across categorical levels. Perform batch analysis.
• Test: Fligner-Killeen homogeneity of variance test: stats::fligner.test(). For when groups lack normal distributions.
• Results: Only Fv/Fm ~ Topography is heteroskedastic (for the photosynthesis Anova’s; ignoring vegetation shade metrics for now).

#batch variance tests for life zone sites (Veg.Type)
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(~ tidy(fligner.test(. ~ D$Veg.Type)),.id="Variable")#put results into a tidy df
kable(V.tests, format = "pandoc", caption = "Site Variance: Fligner-Killeen Equal Variance by Site: T0 summer stress (n=94)", digits= c(0, 3, 4))#plot pretty & round output

#Topography
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(~ tidy(fligner.test(. ~ D$Topography)),.id="Variable")#put results into a tidy df
V.tests<-V.tests %>% select(-parameter, -method) #exclude output parameter & method

#plot pretty
kable(V.tests, format = "pandoc", caption = "Topography Variance: Fligner-Killeen Equal Variance by Topozone: T0 summer stress (n=94)", digits= c(0, 3, 4))#round output

#Microhabitat
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(~ tidy(fligner.test(. ~ D$Habitat)),.id="Variable")#put results into a tidy df
V.tests<-V.tests %>% select(-parameter, -method) #exclude output parameter & method

#plot pretty
kable(V.tests, format = "pandoc", caption = "Habitat Variance: Fligner-Killeen Equal Variance by Microhabitat type: T0 summer stress (n=94)", digits= c(0, 3, 4))#round output

Site Variance: Fligner-Killeen Equal Variance by Site: T0 summer stress (n=94)

Variable	statistic	p.value	parameter	method
Fv.Fm	5.289	0.0711	2	Fligner-Killeen test of homogeneity of variances
PhiPSII	5.627	0.0600	2	Fligner-Killeen test of homogeneity of variances

Topography Variance: Fligner-Killeen Equal Variance by Topozone: T0 summer stress (n=94)

Variable	statistic	p.value
Fv.Fm	8.398	0.0150
PhiPSII	1.271	0.5297

Habitat Variance: Fligner-Killeen Equal Variance by Microhabitat type: T0 summer stress (n=94)

Variable	statistic	p.value
Fv.Fm	1.498	0.4727
PhiPSII	4.100	0.1287

#———————— # ONE-WAY ANOVA’s

• Goal: use a family of univariate 1-way ANOVA’s to test for mean differences in photosynthethetic summer stress (Fv/Fm & Phi-PSII) across categorical predictors: elevation-life zone sites, topography zones (aka topozones), & microhabitat types.
• Unbalanced Design Challenge: missing factor levels (South-slope at Site 1 and interspaces at Site 3) prevent a global 3-way ANOVA. Also, all levels have different sample sizes for all three predictors (another type of unbalanced design).
• Heteroskedasticity Challenge: only Fv/Fm by Topozones is heteroskedastic. For this test I will use Welch’s ANOVA if assumption of normality is met. Otherwise, classic ANOVA is used (whether or not normality is perfect).
• P-values: I explicitly avoid adjusting pvalues for this family of tests because this is exploratory analysis using messy ecological data. We want to report all possible patterns of significance!

Classic ANOVAs

• Tutorial: http://www.sthda.com/english/wiki/one-way-anova-test-in-r#check-the-normality-assumption
• summary.aov() summarizes ANOVA model. Use separately or together in one line. Goal: compute ANOVA & post-hoc multiple comparisons for significant Phi ~ Site ANOVA. Result: Low differs from High Site, but Mid not different from either.
• Effect sizes: using rstatix::anova_test() - produces Cohen’s F, Eta-squared, etc.

#1-way classic ANOVAs - use var.equal = TRUE
W.Aovs<- D %>% select(Fv.Fm, PhiPSII) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Veg.Type, var.equal= TRUE, data= D)), 
                .id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, denom.df
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Classic ANOVAs vs Site - assumes equal variance", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()

Classic ANOVAs vs Site - assumes equal variance

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Fv.Fm	2	91	1.614	0.2048	0.188	0.034
PhiPSII	2	91	2.702	0.0724	0.244	0.056

## Multiple Comparisons
Mod<- aov(PhiPSII ~ Veg.Type, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PhiPSII ~ Veg.Type, data = D, var.equal = TRUE)
## 
## $Veg.Type
##                                   diff          lwr       upr     p adj
## Mid-shrubland-Low-scrubland 0.02457240 -0.058961924 0.1081067 0.7635722
## High-woodland-Low-scrubland 0.07795475 -0.005579571 0.1614891 0.0725895
## High-woodland-Mid-shrubland 0.05338235 -0.024383869 0.1311486 0.2360452

Welch’s ANOVAs: Topography

• use stats::oneway.test() with argument var.equal = FALSE, for Welch Correction, an approximate method of Welch (1951), which generalizes the 2-sample Welch test to the case of arbitrarily many samples. Welch’s correction reduces the denominator degrees of freedom in the F-test.
• Question: I read there are no SS’s in Welch’s ANOVA. Why? *pwr::unknown function but the sjstats package must use it to calculate Cohen’s D!
• Multiple comparisons (Games-Howell): South not different from Flat!

#1-way Welch's ANOVAs
W.Aovs<- D %>% select(Fv.Fm) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= FALSE, data= D)), 
                .id="Variable")#run anova
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract R-sq values from summary.lm(Mod)
R.sq<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-summary.lm(Mod1)
     Out$r.squared
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVA Fv/Fm vs Toopography - heteroskedastic", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()

Welch’s ANOVA Fv/Fm vs Toopography - heteroskedastic

Variable	num.df	den.df	statistic	p.value	Rsq	Cohens_F
Fv.Fm	2	54.745	12.098	0	0.1932	0.489

#Multiple comparisons - Games-Howell
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following objects are masked from 'package:plyr':
## 
##     desc, mutate
## The following object is masked from 'package:stats':
## 
##     filter
games_howell_test(data=D, formula= Fv.Fm ~ Topography, conf.level=0.95, detailed=F)

Residual plots

• Goal: make residual plots of the Welch’s ANOVA model (Fv.Fm ~ Topography) to be fit below to test assumption of normality.
• Results: QQplot looks pretty normal! Assumption satisfied!

#Fv/Fm ~ Topography
Mod1<- aov(Fv.Fm ~ Topography, data=D)
summary.lm(Mod1)
## 
## Call:
## aov(formula = Fv.Fm ~ Topography, data = D)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.288588 -0.081234 -0.001326  0.095766  0.255541 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.43183    0.02823  15.298   <2e-16 ***
## TopographyF -0.04037    0.03595  -1.123   0.2644    
## TopographyN  0.10676    0.03655   2.921   0.0044 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1354 on 91 degrees of freedom
## Multiple R-squared:  0.1932, Adjusted R-squared:  0.1754 
## F-statistic: 10.89 on 2 and 91 DF,  p-value: 5.739e-05
plot(Mod1)

Classic ANOVA: PhiPSII ~ Topography

• use stats::oneway.test() with argument var.equal = TRUE, for classic ANOVA.
• Multiple comparisons result: Phi-PSII differs only between North & Flat!

#1-way Welch's ANOVAs
W.Aovs<- D %>% select(PhiPSII) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= TRUE, data= D)), 
                .id="Variable")#run anova
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract R-sq values from summary.lm(Mod)
R.sq<-D %>% select(PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-summary.lm(Mod1)
     Out$r.squared
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Extract Cohen's F from batch
CF<-D %>% select(PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "ANOVA: Phi vs Toopography - equal variances", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()

ANOVA: Phi vs Toopography - equal variances

Variable	num.df	den.df	statistic	p.value	Rsq	Cohens_F
PhiPSII	2	91	6.483	0.0023	0.1247	0.377

#Multiple comparisons
Mod<- aov(PhiPSII ~ Topography, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PhiPSII ~ Topography, data = D, var.equal = TRUE)
## 
## $Topography
##            diff         lwr        upr     p adj
## F-S -0.06330905 -0.14529288 0.01867479 0.1624954
## N-S  0.04706138 -0.03629741 0.13042018 0.3741169
## N-F  0.11037043  0.03701940 0.18372146 0.0015646

Habitat ANOVAs

• Equal variances, so classic anovas:
• Multiple comparisons (Fv/Fm): Canopy differs from dripline & interspace (interspace = dripline)!
• Multiple comparisons (Phi): All three habitat types differ!!

#1-way classic ANOVAs - use var.equal = TRUE
W.Aovs<- D %>% select(Fv.Fm, PhiPSII) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Habitat, var.equal= TRUE, data= D)), 
                .id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, denom.df
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Classic ANOVAs vs Microhabitat type - equal variance", digits = c(3, 3, 3, 3, 4, 5, 4), table.attr = "style = \"color: black;\"") %>% kable_styling()

Classic ANOVAs vs Microhabitat type - equal variance

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Fv.Fm	2	91	12.433	0	0.523	0.215
PhiPSII	2	91	14.078	0	0.556	0.236

#Multiple comparisons (Fv/Fm)
Mod<- aov(Fv.Fm ~ Habitat, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Fv.Fm ~ Habitat, data = D, var.equal = TRUE)
## 
## $Habitat
##                            diff        lwr         upr     p adj
## Dripline-Canopy     -0.09583631 -0.1715138 -0.02015878 0.0091863
## Interspace-Canopy   -0.17431250 -0.2622704 -0.08635459 0.0000250
## Interspace-Dripline -0.07847619 -0.1746205  0.01766808 0.1321162

#Multiple comparisons (Phi)
Mod<- aov(PhiPSII ~ Habitat, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PhiPSII ~ Habitat, data = D, var.equal = TRUE)
## 
## $Habitat
##                            diff        lwr         upr     p adj
## Dripline-Canopy     -0.11274107 -0.1813238 -0.04415837 0.0005046
## Interspace-Canopy   -0.15511806 -0.2348298 -0.07540627 0.0000349
## Interspace-Dripline -0.04237698 -0.1295077  0.04475369 0.4807655